# Set the working directory
source("~/DALY_function_PSA.R")
library(openxlsx)

#####################
## CIRRHOSIS#########
#####################

#read CSV
# Set the working directory
wd <- paste(wd_root,"/Cirrhosis",sep = "")
setwd(wd)

list_file <- list.files()
list_markov <- list()
name_markov <- c()
for (i in c(1:length(list_file))){
  list_markov[[i]] <- read.csv(list_file[i], header = TRUE, sep = ",", quote = "",row.names("cycle"),
                               dec = ".", fill = TRUE, comment.char = "")
  list_markov[[i]] <- subset(list_markov[[i]],select = c(2:39))
  name_markov[i] <- substr(list_file[i],1,(nchar(list_file[i])-4))
}
names(list_markov) <- as.numeric(gsub("[^0-9]", "", name_markov))

#get cirrhosis incidence
list_Cirr <- list_markov
for (i in c(1:length(list_Cirr))){
  list_Cirr[[i]] <- get_incidence_CC(list_markov[[i]])*1000
}


Cirr <- unlist(list_Cirr, use.names = FALSE)
# distribution

hist(Cirr)
plot(density(Cirr))
#get the mean

Cirr_PSA_mean <- mean(Cirr)
Cirr_PSA_CI95 <- t.test(Cirr)$conf.int
Cirr_PSA_Median <- median(Cirr)
Cirr_PSA_IQR <- quantile(Cirr, probs = c(0.025, 0.975))

Cirr_mean_CI <-  paste(round(Cirr_PSA_mean, 2), " [", round(Cirr_PSA_CI95[1], 2) ,"-",round(Cirr_PSA_CI95[2],2), "]", sep= "")
Cirr_median_IQR <-  paste(round(Cirr_PSA_Median, 2), " [", round(Cirr_PSA_IQR[1], 2) ,"-",round(Cirr_PSA_IQR[2],2), "]", sep= "")

Cirr_recap <- c("Cirrhosis", Cirr_mean_CI, Cirr_median_IQR)
#####################
##    DCC   #########
#####################

#read CSV
# Set the working directory
wd <- paste(wd_root,"/DCC",sep = "")
setwd(wd)

list_file <- list.files()
list_markov_DCC <- list()
name_markov_DCC <- c()
for (i in c(1:length(list_file))){
  list_markov_DCC[[i]] <- read.csv(list_file[i], header = TRUE, sep = ",", quote = "",row.names("cycle"),
                               dec = ".", fill = TRUE, comment.char = "")
  list_markov_DCC[[i]] <- subset(list_markov_DCC[[i]],select = c(2:39))
  name_markov_DCC[i] <- substr(list_file[i],1,(nchar(list_file[i])-4))
}
names(list_markov_DCC) <- name_markov_DCC

#get DCC incidence

list_DCC <- list_markov_DCC
for (i in c(1:length(list_DCC))){
  list_DCC[[i]] <- get_incidence_DCC(list_markov_DCC[[i]])*1000
}

DCC <- unlist(list_DCC)


plot(density(DCC))
DCC_PSA_mean <- mean(DCC)
DCC_PSA_CI95 <- t.test(DCC)$conf.int
DCC_PSA_Median <- median(DCC)
DCC_PSA_IQR <- quantile(DCC, probs = c(0.025, 0.975))

DCC_mean_CI <-  paste(round(DCC_PSA_mean, 2), " [", round(DCC_PSA_CI95[1], 2) ,"-",round(DCC_PSA_CI95[2],2), "]", sep= "")
DCC_median_IQR <-  paste(round(DCC_PSA_Median, 2), " [", round(DCC_PSA_IQR[1], 2) ,"-",round(DCC_PSA_IQR[2],2), "]", sep= "")
DCC_recap <- c("DCC", DCC_mean_CI, DCC_median_IQR)

#####################
## HCC ##############
#####################
# Set the working directory
wd <- paste(wd_root,"/HCC",sep = "")
setwd(wd)

list_file <- list.files()
list_markov_HCC <- list()
name_markov_HCC <- c()
for (i in c(1:length(list_file))){
  list_markov_HCC[[i]] <- read.csv(list_file[i], header = TRUE, sep = ",", quote = "",row.names("cycle"),
                               dec = ".", fill = TRUE, comment.char = "")
  list_markov_HCC[[i]] <- subset(list_markov_HCC[[i]],select = c(2:39))
  name_markov_HCC[i] <- substr(list_file[i],26,(nchar(list_file[i])))
}
names(list_markov_HCC) <- as.numeric(gsub("[^0-9]", "", name_markov_HCC))

#get global HCC incidence
list_HCC <- list_markov_HCC
for (i in c(1:length(list_HCC))){
  list_HCC[[i]] <- get_incidence_HCC(list_markov_HCC[[i]])*1000
}


HCC <- unlist(list_HCC)

plot(density(HCC))
HCC_PSA_mean <- mean(HCC)
HCC_PSA_CI95 <- t.test(HCC)$conf.int
HCC_PSA_Median <- median(HCC)
HCC_PSA_IQR <- quantile(HCC, probs = c(0.025, 0.975))

HCC_mean_CI <-  paste(round(HCC_PSA_mean, 2), " [", round(HCC_PSA_CI95[1], 2) ,"-",round(HCC_PSA_CI95[2],2), "]", sep= "")
HCC_median_IQR <-  paste(round(HCC_PSA_Median, 2), " [", round(HCC_PSA_IQR[1], 2) ,"-",round(HCC_PSA_IQR[2],2), "]", sep= "")
HCC_recap <- c("HCC", HCC_mean_CI, HCC_median_IQR)


recap_dir <- paste(wd_root,"/recap HCC Cirr Inc",sep = "")
dir.create(recap_dir)
Recap <-cbind(c("","mean and 95IC", "median and IQR"),Cirr_recap,DCC_recap,HCC_recap)
write.csv2(Recap,file = paste(recap_dir,"/recap.csv",sep = ""))
